#!/usr/bin/perl -w
#Purpose: Open a file, read down, count peaks (zero values) and windows (all non "-" values), return p
#Edit: Now allow multi gene files
#Edit (1/27/2010): Refactor for understanding and code cleanliness
#Usage: pvaluelooper.pl <infile>

my $numargs = $#ARGV + 1;
if ($numargs != 1) {
	print STDERR "Usage: pvaluelooper.pl <infile>\n";
}
$filename = $ARGV[0];
$fileprefix = ">d-";
open(FILENAME,$filename);
open(OUTFILE,">$filename.out");

sub calccolpvalues {
	while(($inputline = <FILENAME>) && ($inputline !~ /^$/)) {
		$empties = 0;
		chomp $inputline;
		@checkarray = split(/\t/,$inputline);
		for($m=0; $m<$arraylength; ++$m) {
			if ($checkarray[$m] eq "-") {++$empties;} 
			elsif ($checkarray[$m] > 0) {++$countwin[$m];}
			else {++$countwin[$m];++$countpeak[$m];}
		}
	}
	$printer = join("\t",@countpeak);
	$printers = join("\t",@countwin);

	@pvalues = ();
	for($m=0; $m<$arraylength; ++$m) {
		$pvalues[$m] = ($countpeak[$m]/$countwin[$m]);
		}
	$printert = join("\t",@pvalues);
}

sub calccellpvalue {
	my ($curdist, @probs) = @_;
	my @pcalc = ();
	foreach $prob (@probs) {
		$ranged = (2 * $curdist) + 1;
		$value = 1-((1-$prob)**$ranged);
		push(@pcalc, $value);
	}
	return @pcalc;
}

sub calcrowpvalue {
#Takes the following variables:
#@distance - the values of the row
#@probability - the frequency of peaks for each row
	print STDERR "Row is: " . join("\t", @distance) . "\n";
	@distanced = sort {$a <=> $b} @distance;	
	$dimvalue = scalar @distanced;
	$masterprob = 1;
	for($q=0; $q<$dimvalue; ++$q) {
		@pcalc = &calccellpvalue($distanced[$q], @probability);
		$semiprob = &bernoulli_exact($q, @pcalc);
		if($semiprob < $masterprob) {$masterprob = $semiprob;}
	}
	return $masterprob;

}

sub bernoulli_exact {
#Takes the following variables:
#@pcalc - the probabilities for each independent variable
	my ($tgt_num, @pcalc) = @_;
	my $dimvalue = scalar @pcalc;
	$semiprob = 0;
	$maxvalue = 2**$dimvalue;
	@binarylist = ();
	for ($runner = 0; $runner < $dimvalue; ++$runner) {
		push(@binarylist,0);
	}
	$binarylist[0] = -1;
	for($n = 0; $n < $maxvalue; ++$n) {
		$sum = 0;
		++$binarylist[0];
		for ($runner = 0; $runner < $dimvalue; ++$runner) {
			$uppity = $runner + 1;
			if($binarylist[$runner] == 2) {
				$binarylist[$runner]=0;
				++$binarylist[$uppity];
			}
		}

		foreach $individual (@binarylist) {$sum = $sum + $individual;}
		$microprob = 1;
		for($p = 0; $p < $dimvalue; ++$p) {
			if($binarylist[$p] == 0) {$z = 1-$pcalc[$p];} else {$z = $pcalc[$p];}
			$microprob = $microprob * $z;
		}
		#if there are as many or more peaks as in the actual row then we add that probability to the total probability 
		if($sum>$tgt_num) {$semiprob = $semiprob + $microprob;}
	}
	if($semiprob == 0) {$semiprob = 1;}
	return $semiprob;
}


#begin main loop of batch program
while (my $firstline = <FILENAME>) {
	#get genename for the current ortholog set
	$genename = [split(/\t/, $firstline)]->[0];
	print OUTFILE "$genename\n";
	# This first line is the names of the genes
	$firstline = <FILENAME>;
	chomp $firstline;
	@setuparray = split(/\t/,$firstline);
	$arraylength = scalar @setuparray;
	@countpeak = ();
	@countwin = ();
	my $startpos = tell(FILENAME);
	# Create empty arrays of the proper length
	for($m=0; $m<$arraylength; ++$m) {
		$countpeak[$m] = 0;
		$countwin[$m] = 0;
	}
	&calccolpvalues();
	foreach my $pval (@pvalues) {
		print STDERR $pval . "\t";
	}
	print STDERR "\n";
	$firstline = <FILENAME>;
	seek(FILENAME, $startpos, 0);

#RORY - The following segment operates on each row.  It's clear that this portion is independent and therefore readily parallelizable 
	while(($inputline = <FILENAME>) && ($inputline !~ /^$/)) {
		chomp $inputline;
		$empties = 0;
		@checkarray = split(/\t/,$inputline);
		@probability = ();
		@distance = ();
		# clean out columns without data 
		for($m=0; $m<$arraylength; ++$m) {
			if ($checkarray[$m] eq "-") {
				++$empties;
			} else {
				push(@probability,$pvalues[$m]); 
				push(@distance,$checkarray[$m]);
			}
		}
		print OUTFILE &calcrowpvalue();
		print OUTFILE "\n";
	}
}
close (OUTFILE);
close(FILENAME);
exit;
